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Abstract. We propose a numerical approach to study the mechanics of a flowing 
bubble in a constraint micro channel. Using an open source two phase flow solver (Gerris, 
gfs.sourceforge.net) we compute solutions of the bubble dynamics (i.e. shape and terminal 
velocity) induced by the interaction between the bubble movement, the Laplace pressure 
variation, and the lubrication film near the channel wall. Quantitative and qualitative 
results are presented and compared against both theory and experimental data for small 
Capillary numbers. We discuss the technical issues of explicit integration methods on 
small Capillary numbers computations, and the possibility of adding Van der Walls forces 
to give a more precise picture of the Droplet-based microfluidic problem. 


1 INTRODUCTION 

Droplet-based microfluidics is a very promising tool for performing biochemical or 
chemical assays. Droplets are unit systems of controlled volume and content, within which 
mixing can be easily achieved. Several physical phenomena (mechanics, thermocapillarity, 
solutocapillarity, thermomechanics) either in cumulative or compensative ways appears 
when we develop microflnidic setups. It is of prime importance to characterize, under 
controlled experimental conditions, within which range each contribution is the dominant 
phenomena regarding element migration. Rationalizing these various effects would have 
important conseqnences for lab-on-a-chips, and numerical studies are an interesting way 
to understand each contribution separately. 
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In microfluidic setups we often have measurements of shape deformations, bubble or 
drop velocities, and pressures or flow rates at the input/exit conditions, but the knowledge 
local values of these variables are not easily available, because principally of the small 
length scales involved in the system. Numerical approaches are then an interesting way for 
acceding to small length dynamical helds. The validity of numerical approaches requires 
the validation and the confrontation against theories and experimental data. Theoretical 
results exist in few academic microfluidic conhgurations like bubble in cylinders or plates, 
and it is then necessary to be able to compare positively these academic conhgurations 
before extend the prediction to more complex situations which are actually common in 
microhuidic devices. 

In this communication we present an open source two phase how tool (gfs.sourceforge.net) 
for computing solutions of the bubble dynamics (i.e. shape and terminal velocity) in mi¬ 
crohuidic channel and; numerical results are presented and compared, quantitatively and 
qualitatively, against both theory and experimental data for small Capillary numbers. We 
discuss the possibility of adding Van der Walls forces to give a more precise picture of the 
Droplet-based microhuidic problem. 

2 EQUATIONS AND NUMERICAL SCHEME 

We use the incompressible, two-dimensional variable-density, NavierStokes equations 
with surface tension which can be written 

dU 

p{— + UVU) = -Vp + ^iV^U + aK5n ( 1 ) 

VU = 0 

with U = (m, v) the huid velocity, p = p{x, t) the huid density, p = /i(x, t) the dynamic 
viscosity. The Dirac distribution function 5 expresses the fact that the surface tension 
term is concentrated on the interface; a is the surface tension coefficient, k and n the 
curvature and normal to the interface. 

For two-phase hows we introduce the volume fraction c(a:, t) of the hrst huid and dehne 
the density and viscosity as a function of c, i.e. p = p{c{x,t)) and p = p{c{x,t)). The 
advection equation for the density can then be replaced with an equivalent advection 
equation for the volume fraction 


dtC + V{Uc) = 0 

The Navier-Stokes equations are solved using a hnite volume approach based into a 
projection method. The numerical problem is solved using the open-source package Gerris 
[0]. A staggered in time discretisation of the volume-fraction/density and pressure leads 
to the a formally second-order accurate time discretisation. The interface between the 
diherent huids are tracked and followed using a VOF (Volume of Fluid) method. The 
spatial discretisation is done using a quad-tree square cells which give a very important 
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flexibility allowing dynamical costless grid refinement into user-defined regions. Finally 
a powerful discretisation scheme was developed to capture accurately the surface tension 
term. More information can be found in reference [1] 

3 BUBBLE IN A CHANNEL 

The studied configuration is presented in Figure a bubble is pushed into a mi- 
crochannel of width FT by a mean flow velocity Uf. The typical length of the bubble is 
larger than the width H, the bubble is then constrained by the channel. In the stationary 
regime the velocity bubble is Ub- Using the width H as characteristic length and the mean 
flow velocity Uf as characteristic velocity the dimensionless Navier-Stokes is then written 


Uf 



Figure 1: Microfluid configuration : channel width is H, velocities are Uf the mean flow velocity and 
Ub the bubble velocity. 


dU 


+ UVU 


-Vp+Lv^c/+.^ 

lLp_ ILpX^f] 


nSn 


( 2 ) 


where all dynamical variables are dimensionless and we define the Reynolds and Cap¬ 
illary numbers as 


Rf = 


pUH 


a = 


u 

uUf 

a 


(3) 

(4) 


Over this communication we do not discriminate densities and viscosities for liquid and 
bubble as long as in the subsequent computations we £x the density and the viscosity ratios 
between the fluid and the bubble to one, then pf/pb = 1 and Pf/pb = 1- Bretherthon 
[2] studied theoretically and experimentally the dynamic of a bubble on a cylindrical 
conhguration. In the limit of small capillarity number based on the bubble velocity, it is 
shown that the ratio between the gap of the thin film of lubrication h (between the wall 
and the bubble) and the typical height Hof the channel as well as the ratio between the 
bubble and mean fluid velocity scale both as : 


h Ub 

- r\j - 

H Uf 



(5) 
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the proportionality constant depends, at least, on the geometrical conhguration (planar, 
squared or cylindrical) and on the viscosity ratio. By inspecting the h/H relation and 
its dependency into the capillary number Ca, it appears that the grid rehnement plays 
an important role if we want to be able to capture the dynamics of the thin him. The 
Bretherton theory stands that thin him is very important as long as is the key point 
determining the bubble shape. The Figure (left) present the computed shape for a 
capillary number of 0.01, which is in fully agreement with the theory. The Figure 
(right) shows the detail of the grid rehnement at the rear of the bubble where the him is 
thinner. We can also note the specihc grid rehnement along the interface. 
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Figure 2: (left) Typical bubble in the stationary state, (right) Grid refinement at the rear of the bubble. 


The bubble velocity is evaluated by computing the x position of the center mass of the 
bubble along the channel, where x is the spatial position of each fraction c. 

We present now some quantitative and qualitative numerical results on a bubble howing 
on a micro channel. For a bubble between parallel plates the analytical solutions are 


^ ~ 1 + 0.643 (3C'a)2/^ or ~ 1 + 0.51 (OC'a)^/^ (6) 

Uf 


the second relation valid for very viscous drops (liquid-liquid interfaces). In our sim¬ 
ulations we impose a Reynolds number of 0.1 and the only variable parameter is the 
capillarity number Ca- The Figure shows the log-log scaling of the excess of velocity 
^ — 1 as function of the capillarity number Ca up to a capillarity number of 0.5 10“^. 
The solid lines shows both limits from the later relations. These results are consistent 
with those of reference [3] where capillarity number are indeed greater. 

Experimental observations in a microfluidic setup [3] show 

• that below a capillarity number of around 10“® the him gap and the bubble velocity 
remain constant giving no dependency of the ratios h/H and Ub/Uf on the capillarity 
number, and 

• the bubble lost its Bretherton shape becoming more symmetric, like a pancake. 
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Figure 3: Log-log scaling : ^ — 1 as function of Ca- 

This phenomenon appears for h of the order of tens of nanometers, an argument advanced 
as explanation is that in this region, near of the wall, the Van der Waals forces are not 
negligible pushing away the bubble by the apparition of an equilibrium film of constant 
width. The main difficult for including the Wan der Waals forces in this continuum 
approach is that’s necessary resolving numerically a lot of scales, from the small, few 
nanometers for the equilibrium film, to the large ones, the channel width H. Including 
Van der Waals forces in a continuum approach was recently done using Gerris to impose 
the macroscopic contact angles from microscopic physics |Tj, the numerical computations 
were done locally and the wide scale range problem was not matter of fact. Starting from 
the Lennard-Jones potential of two particles and doing some approximations we can add 
to the r.s.h. of the Navier-Stokes momentum equation (equation ([^) the force F{d) per 
unit of volume which depends only on the distance d between the bubble interface and 
the wall 


F{d) 


K 

d* 



(7) 


where it' is a constant and h* is the equilibrium him thickness, m = 3 and n = 2. (details 
in reference H!)- 

To compare qualitatively this approach without resolving all the spatial scales we 
impose a large value of h* which is indeed not physical but the mechanism is still the 
same, pushing the bubble away from the wall. The Figure present two numerical 
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Figure 4: Steady bubble shape and final position including (upped side) or not (lower side) the Van der 
Waals forces. 

simulations of a bubble flowing into a channel for Ca = 0.01 and Re = 0.1 for the same 
hnal time with (upped side) and with out (lower side) Van der Waals forces from equation 
([^. We observe that (i) imposing the Van der Waals forces we found a larger gap h and 
consequently a faster bubble velocity Uh, (ii) the bubble shape becomes a pancake like. 

4 CONCLUSIONS 

We have presented numerical simulations of a bubble into a channel, the well behavior 
of the numerical implementation of the Navier-Stokes equations with a surface tension 
model was demonstrated by a comparison with Bretherton theory for very small capillary 
numbers were the scaling law in was validated. The quality of the numerical results 
are, in particular, a consequence of the grid rehnement approach which allows computing 
the very thin dims of liquid between the bubble and the wall. We have also implemented 
a Van der Waals like force and the imposed minimum gap gives numerical prediction in 
according with experimental observations. 
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